Benchmark of a modified Iterated Perturbation Theory approach on the FCC lattice 

at strong couphng 
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The Dynamical Mean-Field theory (DMFT) approach to the Hubbard model requires a method to 
solve the problem of a quantum impurity in a bath of non-interacting electrons. Iterated Perturbation 
Theory (IPT) has proven its effectiveness as a solver in many cases of interest. Based on general 
principles and on comparisons with an essentially exact Continuous- Time Quantum Monte Carlo 
(CTQMC) solver, here we show that the standard implementation of IPT fails away from half-filling 
when the interaction strength is much larger than the bandwidth. We propose a slight modification 
to the IPT algorithm that replaces one of the equations by the requirement that double occupancy 
calculated with IPT gives the correct value. We call this method IPT-D. We recover the Fermi 
liquid ground state away from half-filling. The Fermi liquid parameters, density of states, chemical 
potential, energy and specific heat on the FCC lattice are calculated with both IPT-D and CTQMC 
as benchmark examples. We also calculated the resistivity and the optical conductivity within IPT- 
D. Particle-hole asymmetry persists even at coupling twice the bandwidth. Several algorithms that 
speed up the calculations are described in appendices. 



I. INTRODUCTION 

Within the last fifteen years or so, the Dynamical 
mean field theory approach (DMFT)^"'^ and its cluster 
generalizations'*"^ have become some of the most power- 
ful techniques to study strongly correlated electrons. In 
these approaches, a single-site hybridized to a bath or a 
cluster hybridized to a bath must be solved. The bath of 
non-interacting electrons is determined self-consistently. 
At the heart of the DMFT approach then, one finds so- 
called impurity solvers. There are now very powerful im- 
purity solvers, for example Continuous Time Quantum 
Monte Carlo (CTQMC) methods^. These methods are 
exact within statistical errors and, for the one-band Hub- 
bard model, certain versions^ at the single site level do 
not suffer from sign problem. Yet, these approaches re- 
quire sizeable computational resources and, in addition, 
real frequency information must be obtained through 
analytical continuation of data with statistical uncertain- 
ties, an ill-posed problem^". It is thus still of great interest 
to work with approximate solvers that are reliable and do 
not suffer from statistical uncertainties. This facilitates 
the calculation of real-frequency quantities with Fade 
approximants^* or directly in real-frequency and also en- 
ables one to quickly explore phase diagrams and pin- 
point interesting regions of parameter space where state 
of the art solvers would be useful. Among possible ap- 
proximate solvers, one finds exact diagonalization, slave 
bosons, Non-Crossing Approximation (NCA), Numeri- 
cal Renormalization Group (NRG) and others^. They all 
have advantages and disadvantages. For example, exact 
diagonalization can consider only a limited number of 
bath sites, NCA is limited to high temperatures and NRG 
to low energies. 

Here we consider Iterated Perturbation Theory 
(IPT)-"^^, an interpolation approach that generalizes the 
original^^^^'* IPT applicable only at half filling. This me- 



thod has been, and still is, wildly used^'^^. The inter- 
polation is constructed so that the self-energy recovers 
both the exact result in the atomic limit and the high- 
frequency limit of the Hubbard model. There is one pa- 
rameter however that cannot be determined from these 
constraints. There have been several proposals to fix this 
parameter. At T = one can impose that Luttinger's 
theorem be satisfied (IPT-L) as was done in [12] but when 
this condition is applied at finite temperature, the results 
are not satisfactory*®. Another very popular approach for 
non-zero temperature fixes the occupation hq of the non- 
interacting part of the Anderson impurity problem used 
in the perturbative calculation to be equal to the lattice 



occupation n 



no 
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(IPT-no)- This condition is ar- 



bitrary since there is no general principle relating these 
two numbers, but it turns out to be quite satisfactory in 
the case of correlated metal i.e. U < Umou^'^'^^'^^ ■ Umou 
is the coupling for which the metal to insulator Mott 
transition occurs at half-filling. 

Despite this success, it is known*^ that when U is larger 
than the critical value for the Mott transition at n = 1 
[U > Umou), then IPT breaks down for n > 1 at low 
T when the condition n = tiq is applied. This happens 
even if, in principle, IPT is constructed to respect the 
atomic limit U ^ t. li has been proposed**" that preser- 
ving the third moment of the spectral weight improves 
the results. Here we show that IPT-np is unsatisfactory 
for U >• Umou close to half-filling for both n > 1 and 
n < 1. We propose a way to circumvent this problem by 
using the fact that when U is large enough, the double 
occupancy becomes almost temperature independent in 
the paramagnetic state with a value that is, to a high 
degree of accuracy, a simple function of the density. This 
provides us with a condition different from n — uq that 
allows one to close the IPT equations even for large cou- 
pling. This approach, IPT-D, is applicable at all tempe- 
ratures contrary to the approach that enforces Luttin- 
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gcr's theorem. It can in principle be improved further by 
enforcing the third- moment sum rule^^. 

In Sec. II we summarize the DMFT approach, the sol- 
vers that we use and the manner in which Fermi liquid 
parameters are extracted. Sec. Ill demonstrates the fai- 
lure of IPT at large coupling. In this section and throu- 
ghout the text, numerical examples are obtained with 
the 3-dimensional FCC lattice. Amongst lattice presen- 
ting electronic frustration, the FCC lattice is important 
because of its prevalence in nature. Our main contribu- 
tion appears in Sec. IV where we show that double oc- 
cupancy can be accurately determined from simple argu- 
ments at very strong coupling and then used to fix the 
remaining parameter in IPT. We call this approach IPT- 
D. Fermi liquid parameters, density of states, chemical 
potential, energy, and specific heat on the FCC lattice 
are calculated with both IPT-D and CTQMC as bench- 
mark examples. Resistivity and optical conductivity ob- 
tained with IPT-f are physically reasonable. Appendix 
A contains details on the three dimensional adaptive inte- 
grator we developed for both IPT and CTQMC calcula- 
tion. Appendix B contains details of the implementation 
of IPT-D. Appendix C explains how to calculate the dif- 
ferent non-interacting functions and Appendix D gives 
details on the calculation of the optical conductivity. 



II. MODEL, DMFT AND IMPURITY SOLVER 
A. Model and DMFT 

We study the one-band Hubbard model, 

H = -^ ti,jd\,^dj,a + UY^ Ui^nn (1) 

where ti,j is the hopping matrix between sites i and j, U 
the on-site Coulomb repulsion, cZ-^J the creation (annihi- 
lation) operator for an electron of spin a on site i and 
nia = d\ ^di^a is the number operator. 

The dynamical mean-field theory (DMFT) provides a 
solution of the Hubbard model that describes the Mott 
transition in three dimensions and has predictive power 
for real materials^'^. Drawing from ideas on the solu- 
tion of the Hubbard model in infinite dimension^^, the 
self-energy in this approach depends only on frequency. 
One first solves the problem of a single site with the 
Hubbard U , hybridized with an infinite bath of non- 
interacting electrons, the so-called Anderson model. One 
extracts the frequency-dependent self-energy of the An- 
derson model, which is then taken as the self-energy in 
the lattice Green's function. The bath is determined self- 
consistently by requiring that projection of the lattice 
Green's function on a single site is identical to the single- 
site Green's function of the Anderson model. The An- 
derson impurity problem can be solved numerically with 
a very high precision. DMFT has been justified with a 
variety of approaches'^ including a variational one^^. The 



single-site DMFT is exact in infinite dimension^. Bench- 
marks against the Bethc ansatz solution in one-dimension 
shows that DMFT can be an accurate solution of the 
Hubbard model also in lower dimensions^^'^^. 

Mathematically, the partition function for the Ander- 
son impurity problem is given by the imaginary-time 
Grassmann path integral 



(2) 



with h = = ksT, A the bath hybridization and 

Sq the action of the impurity, which consists of a single 
site with repulsion U . The self-consistency condition in 

Matsubara frequency reads 



A(iu;„) = iujn + M - S(ia;„) 
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H -1 



- /X - £k - S(iw„) 



(3) 



with E the self-energy. On the FCC lattice the single- 
particle dispersion is given, with lattice spacing a = 
1, by £k = —^t[cos{kx) cos{ky) + cos{kx) cos{kz) + 
cos{ky) cos{kz)]- In single-site DMFT, the self-energy is 
local (I](cj)) which enables one to transform the integrals 
over the Brillouin zone entering the self-consistency re- 
lation into integrals over the non-interacting density of 
states No{oj). However, in the case of a 3d FCC lattice, 
there is no analytic form for Nq (cj) and its accurate nume- 
rical calculation using a Monte- Carlo binning procedure 
is cumbersome and only produces a fixed finite numbers 
of points. On the other hand, if we used a Lorentzian as 
an approximation for the delta function, the band-edges 
and the Van-Hove singularities would suffer from accu- 
racy problems that could be transferred to the DMFT 
calculation. We thus performed the calculation with the 
full k space integration. We devised an adaptive 3d fifth 
order Gaussian quadrature for a cube. This integration 
method is explained in Appendix A. For comparisons 
in calculation of transport properties, we have neverthe- 
less calculated No{oj) using Monte Carlo integration as 
explained in Appendix C. The resulting non-interacting 
density of state for the FCC lattice with nearest-neighbor 
hopping only is shown in Fig. I. 

We have used two "impurity solvers" for the auxiliary 
Anderson model. They are described in the following sub- 
sections. 



B. CTQMC 

The first method is the numerically exact continuous 
time quantum Monte Carlo method (CTQMC)^^, a fi- 
nite temperature approach that relies on the Monte Carlo 
summation of all diagrams obtained from the expansion 
of the partition function in powers of the hybridization A. 
This method does not have a sign problem, and does not 
have errors associated with time discretization or bath 
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Figure 1: Non-interacting density of states for the FCC lat- 
tice with nearest-neighbor hopping only. The large particle- 
hole asymmetry caused by frustration is apparent. 



parametrization. It is therefore exact within statistical 
errors but computationally expensive. We refer to the li- 
terature for an explanation of the approach®'^. 



C. IPT 

We describe the second approach, Iterated Perturba- 
tion Theory (IPT), in more details since it is the focus of 
this paper. IPT is an approximation method that relies 
on an interpolation from 2"*^ order perturbation theory 
for the Anderson impurity problem^^. The interpolation 
preserves the correct high-frequency limit for the self- 
energy and is exact in both the non-interacting and the 
atomic limits. We only consider paramagnetic solutions. 

The self-energy in this approach is parametrized by 



(4) 



where 



E(2)(ia;„) = t/2 f dTe''^-^GZ{T)G^%-T)G^^{T), (5) 
Jo 



with 



Go(«w„) 



1 



iujn + Mo - A(iw„) 



(6) 



and A the hybridization function. The constants A and 
B 



A = 
B = 



n(2 - n) 
no(2 - no) 

(l-f)^ + /io-/i 



(7) 



where no = 26*0 (r = ) and n = 2G(t = ), arc 
chosen such that one recovers the exact solution in the 



atomic limit as well as the exact result for arbitrary U 
in the high-frequency limit. The Green's function used to 
obtain the density n is 



G(iw„) = ^ 



1 



lUlr, 



(8) 



In Eq. (7), ji is the chemical potential of the lattice that 

is determined by fixing the value of n while /io is the 
chemical potential determined by the fictitious density 
no. 

Wc need an additional equation to fix /io. This problem 
has been studied carefully in Refs.[17,19]. Setting fj. — hq 
is not a good option. Indeed, as mentioned in the in- 
troduction, fixing Luttinger's volume works only at very 
low temperature^®. A widely used approach^^ consists in 
fixing n = hq. We call this approach IPT- no. One can 
also modify the formula for the interpolated self-energy 
by requiring that the third moment, appearing in the 
high-frequency expansion of the Green's function, be sa- 
tisfied exactly. In this case, the deficiencies of IPT- no at 
strong coupling are not as severe. We will see below that 
requiring that double-occupancy takes its exact value is 
an easier solution that does not modify the simplicity of 
the original scheme and gives accurate results. 

IPT can be implemented efficiently, as described in Ap- 
pendix. B, so that the solution can be obtained in a very 
short time. 



III. BREAKDOWN OF IPT 

In this section, we first define the physical parameters 
that will be used to demonstrate the breakdown of IPT- 
no. Then wc take advantage of the existence of the exact 
CTQMC impurity solver to characterize the solution of 
the DMFT equation. The last subsection demonstrates 
that for U much larger than the bandwidth, IPT-no fails 
to reproduce even qualitatively the exact solution. 



A. Extracting the Fermi liquid parameters 

At low T and finite doping, it is known that DMFT 

predicts a Fermi liquid regime no matter how strong the 
interaction^. In other words, a quasiparticle peak always 
appears at w = at low T except at half-filling when 
U > Umou- We characterize the Fermi liquid with three 
parameters namely the effective chemical potential. 



/i = /x-I]'(0), 



the quasiparticle weight, 



5I]'(w) 



dui 



(9) 



(10) 



and the scattering rate S"(0), where real and imaginary 
parts are defined by S = E' iS". In the DMFT treat- 



4 



ment of the Hubbard model, Luttinger's theorem is sa- 
tisfied at T = when ft takes the value of the non- 
interacting chemical potential that gives the same den- 
sity. 

All of the above parameters are calculated with the 
self-energy on the real frequency axis and thus, in prin- 
ciple, one needs to perform an analytical continuation 
from the data in Matsubara frequencies and then extra- 
polate to zero temperature. In practice, we calculate the 
values of the self-energy for a few very low temperatures 
and use them to extrapolate to zero frequency and zero 
temperature. For the retarded S'(0) and S"(0) we take 
S(w„=o) at the smallest positive uJn=o for three low tem- 
peratures and extrapolate to T = using the fact that 
i;(w„ 0+) = 0). For Z, the spectral defini- 

tions of the self-energies 



E'(c.) 
S(«w„) 



du' S"(a;') 



TT UJ 

duj' E" 



(11) 



TT UJ' 



allow one to prove 
Im[I](ia;„)] 



da;'E(w') dYI [uj) 



UJ'- 



duj 



{12) 



also 



which, for a linear dependence of Im[I](za;„)] on ujn 
follows from the Cauchy-Riemann relation for holomor- 
phic functions of a complex variable. This last equa- 
tion with three low temperatures allows us to calculate 
Im[I](ia;„=o)]/'^rt=0; extrapolate to T = and obtain Z. 



B. Expected behavior, as obtained from CTQMC 

Consider the one-band Hubbard model on the FCC lat- 
tice, where the single particle dispersion is given by = 
—^t[cos{kx) coii{ky) -f cos(fca;) cos(fcz) -t- cos(fcj^) cos(fc^)]. 
We present the DMFT results obtained with the CT- 
QMC impurity solver for U — 8t and U — 32t, below and 
above Umou for the Mott transition at half-filling. The 
bandwidth is 16i for the 3d FCC lattice. 

Fig. 2 displays the Fermi liquid parameters. In Fig.2- 
(a) the red solid line shows the non-interacting chemical 
potential as a function of density. The effective chemical 
potential jl = /i — E'(0) is shown with blue circles for 
U = 8t and black points for U = 32t. The dashed lines 
indicates the position of the band edges for the 3d FCC 
lattice. As expected, except at half-filling for U = 32t, 
Luttinger's theorem is satisfied. In Fig.2-(b), the extra- 
polated scattering rate E"(a;) is negligibly small, except 
for U = 32t at n = 1. For U — 8t, (not shown) it is of 
the order 1x10^'* and has essentially no density depen- 
dance. The value of Z, shown in In Fig.2-(c), behaves 
as expected : For U > Umou, Z vanishes when the oc- 
cupation approaches half-filling while it is close to the 
non-interacting value Z = 1 when the lattice is almost 
empty or full. We can also see that even for coupling as 
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Figure 2: (Color online)Results obtained with CTQMC as 
impurity solver are plotted as a function of density and shown 
in blue with circles and line for U = 8 and in black with dots 
and line for U = 32t. In all numerical results, energy units 
are such that t = 1. Boltzmann's constant and the lattice 
spacing are also taken as unity. We obtain the zero-frequency 
limit from a poor man's approach : we take pt = 25, 50 and 
75 and use the value of the function at the lowest Matsubara 
frequency in the three cases to perform the extrapolation, (a) 
Check for Luttinger's theorem : The effective chemical poten- 
tial /i = ^ — E'(0) is equal to the non-interacting chemical po- 
tential shown in red except at half-filling where there is a Mott 
gap for U = 32f. (b) KX, U = 32t the imaginary part of the 
self-energy at zero frequency E"(0) should be zero away from 
half-filling and infinite at half-filling, (c) The single-particle 
spectral weight Z vanishes only a.t n = 1, U = 32t where there 
is a Mott gap. 
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large as U = 32t, the absence of particle-hole symmetry 
in the dispersion relation still leads to a value of Z that is 
not symmetric with respect to half-filling. Clearly, elec- 
tronic frustration plays an important role in the doped 
Mott insulator. 



C. Breakdown of IPT-no 

The IPX equations Eqs. (4)- (8) do not determine the 
value of hq. As mentioned previously, for T = the re- 
quirement that Luttinger's theorem be satisfied (IPT-L) 
provides an additional independent equation, except at 
half-filling for U > Umou- However, Luttinger's theo- 
rem is in general not satisfied at finite temperature and 
the method becomes inaccurate. The condition n = uq 
has thus been proposed"^'''-^^ (IPT-tiq). It gives satisfac- 
tory results for correlated metals (not for the insulator 
at half-filling). 

The results for the low-temperature extrapolations of 
p, = /i — E'(0) and Z ior U — 8t and U — 32t are shown 
in Fig. 3. Below the Mott transition, U — 8t, the n — uq 
results (brown (*)) are shown. One can detect only a 
very small difference with the solid red line. Luttinger's 
theorem is thus essentially satisfied. 

On the other hand, IPT-no for U = 32t (kaki (□)) gives 
non-physical results not only^'' for n > 1 but, quite gene- 
rally, close to half-filling. Not only is Luttinger's theorem 
strongly violated, but for a large range of densities, n > 1, 
jl is outside the band. Many properties of the Fermi li- 
quid are proportional to functions of the non-interacting 
system evaluated at fl. But these functions are zero out- 
side the band and so if jl is outside the band we obtain 
zero. For example, this would predict an insulator away 
from half-filling. The situation is not better for Z, es- 
pecially around half-filling where it vanishes for a finite 
range of densities when n > 1. This demonstrates that at 
low r, when U is large, IPT-no in its simplest form can- 
not be applied. We must thus search for a new condition 
to explore this region of parameter space. 

IV. IPT DOUBLE OCCUPANCY : IPT-L> 

Imposing exact results such as sum-rules, whenever 
possible, is desirable for any physical theory. Whereas 
the condition n = no is not required by any fundamental 
principle, the self-energy must always obey 

D = |^e^""°^S(*c.„)G(zw„). (13) 

n 

where D =< n-fUi > is double occupancy. Enforcing this 
consistency condition between single-particle properties, 
such as S and G, and a two-particle property, D, has been 
successful in other approaches, such as the Two-Particle- 
Self-Consistent theory.^^'^^. In the regime of interest here, 
strong coupling, D can be accurately estimated and is 
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Figure 3: (Color online) Fermi liquid parameters as a func- 
tion of density. Zero frequency results are obtained with the 
same extrapolation method as in Fig. 2. (a) Check of Luttin- 
ger's theorem. The effective chemical potential fl = fi — 'E'{0) 
should equal the non-interacting value, shown in red, when 
the theorem is satisfied. For U — 8t, the brown asterisks (*) 
obtained with IPT n = no satisfy the theorem. For U — 32t 
results for three different methods are shown : in kaki (□) for 
IPT n = no, in cyan (★) for IPT D„aive and in magenta (0) 
for IPT {D)cTQMC- (b) E"(0) is plotted for U = 32t in ma- 
genta (0) for IPT {D)cTQMC as above, and compared with 
the CTQMC results shown previously in Fig. 2 (black dots joi- 
ned by a line), (c) Quasiparticle spectral weight Z computed 
for different methods and displayed with the same symbols 
as in (a). We compare with the CTQMC results of Fig. 2, 
namely blue symbols (o) with line for U = 8t and black sym- 
bols (.) with line for U — 32t. The results for IPT n = no at 
U — 32t are un-physical since they predict an insulator away 
from half-filling. 
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Figure 4: (Color online) CTQMC results at [/ = 32t for 
double occupancy D — D„aive plotted as a function of tem- 
perature. We define Dnaive ~ for fillings n < 1 and 
Dnaive =71 — 1 for 71 > 1. The various densities are repre- 
sented by different symbols : n — 0.2 (black (o)), 0.4 (blue 
(x)), 0.6 (red (□)), 0.8 (green {<))), 1.0 (yellow (+)), 1.2 (cyan 
(V)), 1.4 (magenta (a)), 1.6 (brown (<)) and 1.8 (kaki (*)). 
The largest deviations from the naive value, occurring close 
to n = 1, are less than 10~^ in absolute value. 
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Figure 5: (Color online) Double occupancy 73 as a function 
of density obtained from CTQMC for U = 32t for three tem- 
peratures : P = 25/t (black (□)), P = 10 /t (blue (o)) and 
P — 0.5/t (red (•)). On this scale, the naive value of D is very 
accurate. The inset is a zoom for densities 7i < 1. 



only very weakly dependent on temperature, as discussed 
in the following subsection. There, we assess the accuracy 
of the approach. 



A. Exact and naive values of double occupancy at 
strong coupling 

For very large U, it is easy to guess that D should de- 
pend only very weakly on temperature. In addition, the 
value of D can be estimated quite accurately. Indeed if 



n < 1 there are unoccupied sites in the lattice and since 
U is large, D should naively be zero. For tt. > 1, D is ne- 
cessarily non-zero. If we start from the half-filled, Mott 
insulator, with one electron per site and add electrons, 
they must go to a site which is already occupied. Thus D 
should simply equal the excess number of electrons mea- 
sured from half- filling, D = tt. — 1. These estimates are 
called Dnaive ■ In reality for n < 1 there are corrections of 
order t/t/ to double occupancy, giving rise to exchange 
interaction, and D is slightly larger than zero. For analo- 
gous reasons, _D for 7i > 1 will always be a bit larger than 
71 — 1 i.e. D ^ {n — 1) + SD. This is obvious for models 
with particle-hole symmetry, but it will be true as well, 
even in the absence of particle-hole symmetry. 

We can verify our estimates with the CTQMC results 
for U = 32t. Within CTQMC, D is calculated directly 
on the impurity by the Monte-Carlo sampling. Fig. 4 dis- 
plays D— Dnaive as a function of temperature for different 
densities. We see that D — Dnaive is small and that the 
T dependence is on the third significant digits. So, for all 
practical purposes, we can assume D to be independent 
of T although it differs from Dnaive ■ The values of D ob- 
tained are very close to the naive expectation but always 
slightly larger. In Fig. 5 we show the double occupancy D 
as a function of the density for three different tempera- 
ture from T = 0.04t to T = 2t. This figure confirms again 
that even if we have some dependence on T, it is quite 
small and the result is a fairly simple function of the den- 
sity. Very similar results have been obtained in Ref. [28] 
in the case of a 3d simple cubic lattice. We thus have a 
rather simple constraint that we can take into account in 
IPT to fix all parameters. We call this approach IPT-D. 
In the next subsection, we will assess the accuracy of this 
approach and verify how the results are modified when 
the exact value of D is used instead of the naive one. 



B. Accuracy of Fermi liquid parameters IPT-D 

We first set D to its naive values, i.e zero for n < 1 
and n—1 for 71- > 1. The results for fl and Z are shown in 
Fig. 3 as cyan stars. Clearly, the value obtained for fi is 
in much better agreement with Luttinger's theorem than 
in the case IPT-uq, although it is still incorrect for densi- 
ties near half-filling. The biggest improvement is that we 
avoid values of fl outside the band. Furthermore, in the 
case of Z, IPT-D is close to the CTMQC results while for 
IPT-no it is quite far^^ from the correct result, leading in 
particular to an un-physical insulator over a finite range 
of densities for 7i > 1. 

As we now show, one can improve the results further 
by using an accurate value of D. That value can be obtai- 
ned from a number of methods, in particular CTQMC. It 
can be computed quite accurately and does not require a 
large number of Matsubara frequencies. Even if CTQMC 
is available, it may be desirable to use IPT-D because the 
calculation can either be done directly in real frequency 
or analytically continued from Matsubara frequencies 
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using simple methods such as Pade approximants^^, whe- 
reas with CTQMC, Maximum entropy^ is necessary. In 
addition, since D has neghgible temperature dependence 
in strong couphng, only one value of D may be sufficient 
with IPT-I? to compute other quantities for a wide range 
of temperatures. 

Since D is not completely T independent, we use an 
average called {D)ctqmc calculated between /? = 75/t 
and /? — 0.5 /t for the purpose of comparison with the 
naive approach. It is calculated from the arithmetic mean 
of the numerical values of D{T). Note that the values of 
D{T) for each n are taken as the arithmetic mean of 
the last four DMFT iterations. The results are shown in 
Fig. 3 as magenta lozenges. We see that the results for jx 
are in excellent agreement with Luttinger's theorem ex- 
cept very close to half- filling where it deviates, but not 
too much. Surprisingly, Z is not as accurate as that obtai- 
ned from the naive estimate of D. But if we look at the re- 
sults for densities between 0.8 and 1.2, the difference bet- 
ween ■^(_D)cTQj\/c ^^'^ ZcTQMC is small and constant and 
^{D) CTQMC correctly extrapolates to zero at half- filling. 
A small difference in D can have a quantitative impact 
on Z, without affecting qualitative trends. For example 
for n = 0.84, the D given by CTQMC is D = 0.00619 ins- 
tead of the naive D = 0, whereas for n = 1.16, CTQMC 
gives D = 0.1637 instead oi D = 0.16. We could imagine 
that because these quantities are at low T, it would be 
better to take a D that is in the low temperature range. 
If we do this, we obtain a {D)ctqmc slightly larger. We 
then find that jl is not really affected while Z is a little 
bit worse than that obtained from the average D over 
the larger T range. Some tuning of D would allow us to 
get a best possible set of jl and Z, but that is clearly 
not the purpose of the exercise. Finally, for E"(0) Fig. 3- 
(b) shows that IPT-Z3 is qualitatively correct while being 
always smaller than CTQMC. 
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Figure 6: (Color online) Density of states for n — 0.84, j3t = 
25 and U = 32t obtained with three methods : black (solid 
line) with CTQMC maxent, blue (- -) with IPT-(D)ctqmc, 
and red (-.) with IPT-Dnaiue. (b) is a zoom of (a) around 
oj = 0. The value at zero frequency is improved when a more 
accurate value of D is used in IPT. 



C. Accuracy of the Density of States and Chemical 
Potential 

It is instructive to look at the Density of States 
obtained from IPT with fixed D and Pade analytical 
continuation^^. We show n = 0.84 and /3 = 25/t as typi- 
cal values in Fig. 6. We compare to the CTQMC values 
obtained from Maximum Entropy analytical continuation 
of G{iujn)- That Green's function is an average over seve- 
ral converged DMFT iterations. The Maximum Entropy 
implementation that was used here is somewhat crude 
and thus we must not really focus on the details. The 
CTQMC results have different errors at different scale, 
i.e very precise at low w„, fiuctuating at intermediate uin 
while at large aj„ the results are analytical. Hence, we 
choose the weight of the entropy term based on heuristic 
considerations, depending on the real- frequency range we 
are interested in. 

As was noted previously^^'^^, in IPT there are states 
in the Mott gap at finite frequency, but their weight is 



small compared to the states everywhere else, namely 
near zero frequency and in the lower and upper Hubbard 
band. Overall, IPT with fixed D compares well with CT- 
QMC, but, at low temperature, what really matters is 
the region near a; = 0. We thus zoom on this region in 
Fig. 6-(b). There, we see that in the vicinity of a; = 0, 
when D = {D)ctqmc, we are quite close to the CTQMC 
values. When we use the naive D, the quasi-particle peak 
is shifted a little bit to the right and so is this why the 
low T results are different even if the shape and values 
of the peaks are similar. 

We also compare the results for an integrated quan- 
tity, /i(T), that is obtained in general by solving n — 
2 J f (lo) p{u))dio with f{uj) the Fermi function and p{uj) 
the density of states. In our case, this is a byproduct of 
the DMFT calculation. No analytical continuation is in- 
volved. In Fig. 7 we show the results for three densities. 
In Fig. 7(a) and (b) we note that the numerical values 
obtained with the different methods differ by at most 
about 10%. The best results are for IPT- (_D)ctqa/c since 
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Figure 7: (Color online) Chemical potential as a function 
of temperature for different IPT approximations, compared 
with the reference CTQMC calculations as black dots with 
line obtained for U = 32t and different densities : (a) n — 
0.80, (b) n = 1.2 and (c) n = 0.84. The three different IPT 
approximations are given in kaki (□) for IPT-no, in cyan (★) 
for IPT-Dnaive, aud in magenta (0) for IPT-(D)ctqmc ■ The 
latter approximation in magenta (0) is best, having a more 
or less doping and temperature independent offset 5jj,/t ~ 0.5 
when compared with the reference CTQMC in black. 



the curves are qualitatively very similar to the CTQMC 
ones with a derivative quantitatively quite close for all 
T. The absolute difference between 1PT-{D)ctqmc and 
CTQMC is almost doping and temperature independent. 
The derivative with respect to temperature for both IPT- 
riQ and IPT-Dnaive is not as good. At high enough tem- 
perature, all methods give similar results. Fig. 7(c) shows 
the result closer to half-filling, comparing CTQMC with 
IPT-{D)cTQMC, the best IPT method. As already dis- 
cussed, IFT-riQ gives un-physical results in the vicinity of 
half-filling. 



D. Energy and Specific Heat 

In this section we compare internal energy and spe- 
cific heat in IPT-D with CTQMC. Within CTQMC, 
energy can be calculated quite accurately with a reaso- 
nable number of Matsubara frequencies. Indeed, it was 
shown by Haule^^ that the kinetic energy (K) is propor- 
tional to the average perturbation order (fc) for a given 
set of parameters. As already discussed, the double oc- 
cupancy is calculated directly by CTQMC and thus the 
total energy is given by 



CTQMC 



(T) = -T(fc) + UD. 



(14) 



In general, for a many-body system, the energy is given 
by the thermal average, in the grand-canonical ensemble, 
of the Hamiltonian. For the Hubbard model we may write 



k.(7 



i 



1 1 



dr 



T=0- 



(15) 



For IPT-_D, we use directly the imaginary time expres- 
sion. Since we generally compute Green's functions only 
for positive imaginary time, we use the equivalent expres- 
sion 



E{T) = + C/n - /i + 1 ^ ekGkir = 0+) 



-E 



9Gfc(T) I 

9t lr=0+' 



(16) 



where we used that for the FCC lattice Efe = 0. Once 
the energy is calculated, the specific heat C„ is given by 
n — rf-E(^) 



dT 
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Figure 8: (Color online) Energy as a function of temperature 
obtained from IPT-D (solid lines) and CTQMC (dashed lines) 
for U = 324. (a) Densities equal to, or below half-filling n = 
0.2 (black (o)), 0.4 (blue (x)), 0.6 (red (□)), 0.8 (green (0)), 
0.84 (cyan (v)), 0.88 (magenta (A)), 0.92 (brown (<)) and 
1.0 (kaki (★)). For densities above half-filling, displayed in (b), 
n = 1.08 (black (o)), 1.2 (blue (x)), 1.4 (red (□)), 1.6 (green 
(0)), 1-8 (cyan (V)), the quantity UDnaive is subtracted from 
the energy to allow the results to fit on the same scale, (c) is 
a zoom for n — 1.08 and (d) a zoom for n = 1.4. 



0.5 1 1.5 2 

T/t 

Figure 9: (Color online) Specific heat at constant filling as a 
function of temperature for U — 32t. (a) Results from IPT-_D 
(solid line) for densities below half-filling : n = 0.2 (black (o)), 
0.4 (blue (x)), 0.6 (red (□)), 0.8 (green (0)), 0.84(cyan (v)), 
0.88 (magenta (a)), 0.92 (brown (<)) and 1.0 (kaki (*)), (b) 
Specific heat from IPT-D (solid line) for densities above half- 
filling n = 1.08 (black (o)), 1.2 (blue (x)), 1.4 (red (□)), 1.6 
(green (0)), 1-8 (cyan (V)), (c) Comparison between IPT-D 
as solid lines and CTQMC as dashed lines for n = 0.6 (red 
(□)), 0.8 (green (0)) and 0.84(cyan (V)). The peak positions 
and shape coincide, even though the absolute values differ. 
In this case, since CTQMC values comme from differentia- 
tion of Monte Carlo data, there is a rather large uncertainty, 
especially for peaks. 
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To obtain values for all temperatures, we have consi- 
dered the fixed D as the one given by the average over 
0.5 < /3 < 75/i of the CTQMC resuhs. This is what 
we called IPT-{D)ctqmc previously. We use IPT-D for 
brevity. We present the comparison in Fig. 8 for different 
densities. Note the change in color code (see the legend). 
In Fig.8-(a) we display results for densities n < 1 while 
in Fig.8-(b), the density is above half-filling n > 1. The 
energy scale for the latter case is shifted by a filling de- 
pendent quantity U Dnaive so that the various fillings can 
be displayed on the same scale. In Fig.8-(c) and (d) we 
zoom on n = 1.08 and n — lA respectively without 
energy shift to emphasize the differences between IPT-D 
and CTQMC. Comparing those differences in Fig.8-(c) 
with those in Fig.8-(d), we see that the further we are 
from half-filling, the better the agreement. 

In Fig. 9 we plot the specific heat as predicted by IPT- 
D. Once again we consider separately n < 1 in Fig.9-(a) 
and n > 1 in Fig.9-(b). The closer we are to half-filling, 
the lower the temperature at which the peak that signals 
the appearance of the Fermi liquid regime appears. Des- 
pite the strong coupling, the particle-hole asymmetry is 
noticeable. 

In Fig.9-(c) we can compare the IPT-Z? and the CT- 
QMC results for three densities n = 0.60, 0.80 and 0.84. 
The CTQMC results for the energy are not completely 
smooth especially at low T because of statistical errors in 
the data. Thus, if we were to calculate the specific heat 
Cn for CTQMC, we would, for many of the values of n, 
have to be very careful before calculating the derivative. 
Here, we just performed a simple derivative on the raw 
CTQMC data to verify the trend. From the results for the 
energy. Fig. 8- (a) we already knew that even if the change 
of curvature in the coherent-incoherent region exists also 
for CTQMC, it is much smoother than the one obtained 
from IPT. This is apparent indeed in C„. Nevertheless, 
the transition temperature for the coherent-incoherent 
crossover given by the position of the maximum is quite 
similar for both methods even if, for CTQMC, its exact 
value is hard to really pinpoint due to the statistical er- 
rors in the raw data. 

We end this subsection by noting that the inadequacy 
of IPT-no also shows when the energy is calculated. One 
finds an increase of the energy with decreasing T at low 
T, which is of course non physical since this corresponds 
to a negative specific heat. 



E. DC resistivity 

Analytical continuation of response functions obtai- 
ned with CTQMC is in general very difficuh. With IPT, 
one can calculate response functions by first analytically 
continuing the self-energy with Pad? approximants and 
using the real-frequency expressions in terms of spectral 
weight. In this section, we obtain results for the elec- 
trical resistivity p within IPT-Z? to verify whether they 
are physically sensible at strong coupling. With IPT-np 
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Figure 10: Transport function X(e) = {w^Y ^(^ - ^fe). 
as calculated in appendix C. 
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Figure 11: (Color online) Resistivity as a function of tempe- 
rature for U = 32t as calculated by IPT-D for diflerent values 
of density : n = 0.2 (black (o)), 0.4 (blue (x)), 0.6 (red (□)), 
0.8 (green (0)), 1.2 (cyan (v)), 1.4 (magenta (a)), 1.6 (brown 
(<)) and 1.8 (kaki (★)). The resistivity are largest close to 
half-filling where they exhibit low coherence temperatures. 



they are not : One obtains an insulator at finite filling on 
the electron-doped side. Since vertex corrections vanish 
in single-site DMFT, the conductivity can be obtained 
from 

(_f!M)^(||)\.„„,, 

(17) 

where <Jo = ^ and A{k,uj) = -ilm {^^^^j^^^i^p^^}. 

We have restored the lattice spacing a to exhibit the 
units. Because the self-energy is local, we can also re- 
place the triple sum over fc by a one dimensional integral 
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over an energy variable. 



cTxxiO) = c^o / deX{e) duj [ - 



du! 



A^{e,Lu), (18) 



where the so called transport function that includes the 

effect of the lattice is X{e) = J2k ^(^ - ^k)- For 

a simple cubic lattice in any dimension, the calculation 
of this function can be brought in the form of a one- 
dimensional integral'^'^ , but on the FCC lattice this is not 
possible. If one wishes to use the expression Eq. (18) ins- 
tead of keeping the full three dimensional integrals in 
Eq. (17), X must be calculated numerically. We explain 
in Appendix C an efhcient way to do this. The result for 
the FCC lattice is shown in Fig. 10. We have tested both 
ways of obtaining axx{0), i-e. Eq. (18) and Eq. (17) and 
they give essentially the same answer. They cannot give 
exactly the same number because X{e) is calculated for a 
fixed number of points and we must thus interpolate bet- 
ween these points when performing the integral over e. 
This adds another source of numerical error not present 
in Eq. (17). But, contrary to the non-interacting density 
of states (Fig. 1), X{e) is a smooth function (Fig. 10) 
and thus the process of interpolation will only induce 
a negligible error. When the integrals are performed in 
the order shown in Eq. (18), i.e. integrate over w first, 
the resulting integrand for the e integration is smooth 
and X(e) need not be obtained with extreme accuracy. 
For these reasons, for the conductivity we preferred to 
use Eq. (18). Using Eq. (17) increases dramatically the 
calculation time here contrary to G(zw„) or xii(*^ra) dis- 
cussed in the next section. 

At low temperature and finite doping, the conductivity 
is proportional to X{il). Thus, if /i is such that its value is 
outside the non- interacting band — 12 . . . 4, X(/i) =0 and 
thus axx{(^) = 0. In an exact implementation of single- 
site DMFT, a null conductivity can only happen at half- 
filling for U > Umou- At any finite doping, there is a 
quasi-particle peak and Luttinger's theorem is respected 
(Fig.2-(a)). As IFT-np fails with respect to Luttinger's 
theorem (Fig.3-(a)) in that case p start diverging at low 
T for the densities with jl outside the band while it should 
exhibit the behavior of a Fermi liquid. 

We show in Fig. 11 the result for IPT-D. We see that 
it has the correct behavior at low temperature. For 
clarity, we omitted values of densities between n = 0.80 
and n = 1.2. Nothing very different happens there. We 
still have the low temperature behavior and, as is 
aheady obvious from the figure, the absolute values of p 
obtained are larger and larger when we approach half- 
filling, the Mott insulating state. Note that close to half- 
filling, even though there is a tendency for the resistivity 
to saturate at very high temperature'^^, this occurs at 
values of the resistivity much larger than the Mott-Ioffe- 
Regel limit p finj . This is characteristic of incoherent 
transport in strongly correlated systems. 
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Figure 12: (Color online) Optical conductivity for U = Zlt 
and n = 0.80, as calculated from IPT-_D for three different 
temperatures using two analytical continuation approaches 
for each temperature. The maximum entropy results are re- 
presented with solid lines and the Pade analytical continua- 
tions with dashed lines. The broadest zero-frequency peak 
(red (x)) is for the largest temperature, /3 = Xjt, and the nar- 
rowest one (black (o)) for the lowest temperature, /3 = 25/t. 
The intermediate case, /3 = 2.?i/t is in blue (□). The features 
in the optical conductivity can be identified with transitions 
between the Fermi level and peaks in the single-particle den- 
sity of states. 



F. Optical conductivity 

We analyzed the performance of our new solver IPT- 
D for Fermi surface properties at T = (Z,I]',I]"), for 
integrated quantities like chemical potential /i, energy, 
specific heat, resistivity and frequency-dependent func- 
tions such as the density of states. To finish, we look at 
the optical conductivity. Appendix. D explains how it is 
calculated using the susceptibility Xii(*^n) in bosonic 
Matsubara frequency and analytical continuation. 

As a first check, which does not depend on analyti- 
cal continuation, we verify the /-sum rule xii(ir2„ = 

0) = Yl,k ^^("fe) (^^® example Ref. [32]). In the case 
of nearest-neighbor hopping, this quantity can be rela- 
ted to the kinetic energy. For a FCC lattice one finds 
Sfc ^^('^fc) = ^|(-^) where {K) is the average kinetic 
energy. For CTQMC the average kinetic energy is obtai- 
ned from {K) = —T{k) where, as already mentioned, (fc) 
is the average perturbation order obtained directly from 
the Monte-Carlo simulation. To calculate xiii^^n — 0) 
one needs the dressed Green's function, or equivalently 
the self-energy (E(ia;„)). A very large number of numeri- 
cal operations is necessary, as explained in Appendix. D, 
but analytical continuation is unnecessary. For all the 



tests we did, the ratio 



xii(in„^o) 



agreed with 2/3 up 



to the third digit. For example, for n = 0.84, P/t = 25 
and U = 32t the ratio is 0.6668, while for n = 0.80 it is 
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0.6664. 

Analytical continuation of CTQMC is problematic, es- 
pecially since we have a very wide frequency range given 
the large value of U. Hence we display only results ob- 
tained with analytical continuation of IPT-I? and check 
for consistency with what is expected from the den- 
sity of states. Results for density n — 0.80, interaction 
[/ = 32 and three temperatures is shown in Fig. 12. 
Solid and dashed lines correspond respectively to Maxi- 
mum Entropy"^^ and Pade analytical continuation. There 
is some quantitative disagreement but the qualitative in- 
formation is the same. At low temperature {^/t — 25 
black lines) there is a clear peak around w « 5 which cor- 
responds to transitions between the lower Hubbard band 
and the quasi-particle peak appearing in the density of 
state for a similar density in Fig. 6. At larger tempera- 
ture (blue lines (3/t = 2.3), the decrease of the peak near 
a; 5 in cr(a;) corresponds to the disappearance of the 
quasi-particle peak and loss of coherence. That loss of co- 
herence for n = 0.80 is signaled by the maximum in C„ 
observed in Fig. 9. At an even larger temperature, (red 
lines 13 /t = 1) the w 5 peak has disappeared since we 
are now in the incoherent regime. The peak for transition 
to the upper Hubbard band around 32i is always visible. 



V. CONCLUSION 

In addition to being numerically inexpensive, IPT pro- 
vides a method where analytically continued results can 
be reliably obtained directly in real frequencies or from 
Pade approximants^^ instead of Maximum Entropy me- 
thods required when Quantum Monte Carlo is used as an 
impurity solver. 

However, for large interaction strengths in doped Mott 
insulators, the popular condition for IPT where one im- 
poses n — no fails at low temperature for a broader re- 
gime than previously expected^^. As a solution, we pro- 
pose that one should instead enforce the exact relation 
Eq.(13) between double occupancy and single particle 
quantities. Further improvements are expected if one also 
enforces the third moment of the spectral weight^^. 

Our new method, IPT-D, can be used for any coupling 
if D is known. Double occupancy D can be obtained quite 
accurately by a number of methods and is negligibly de- 
pendent on temperature for large coupling. For example, 
in the strong coupling regime one can use exact diagona- 
lization of small clusters or slave bosons^ while, at weak 
coupling, methods such as Two-Particle-Self-Consistent 
theory (TPSC)^^'^'' give good results. Also, in both re- 
gimes. Quantum Monte-Carlo methods can be used. In 
the very large U limit, the naive estimate 13 = for 
n < 1 and D — n — 1 for n > 1 leads to qualitatively 
correct results, except very close to half-filling. Our ap- 
proach has been benchmarked on the FCC lattice for a 
number of observables. The large-particle hole asymme- 
try of that lattice survives in observable quantities even 
for interaction strength equal to twice the bandwidth. 
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Figure 13: The special points for a 3d Gaussian quadrature 
of fifth order over a cube of length 2h. 
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Annexe A: Integrator 

As explained in the main text, for accuracy in the 
DMFT iteration we need to perform integrals over a three 
dimensional Brillouin zone. In this appendix we use sym- 
metry and ideas from Gaussian quadrature, adaptive me- 
thods, and statistics, to devise an accurate and fast in- 
tegrator. We first obtain a quadrature of order five and 
then explain how we can make it adaptive. 

We need a normalized triple integral over a cube of 
length 2h centered at a point Tq = [xQ,yQ,ZQ]. If we put 
the origin at this point, the integral takes the form 

1 /"^ 

J^W J I J I J ^^^'^2/c^^-^(a;,y,2^)- (Al) 

We considered a normalized integral because the inte- 
grals we need to solve are over k-space and thus need 
normalization. 
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We first show that fourteen appropriately chosen 
points and only two weights can give us an approxima- 
tion of order five. Usually when one develops a Gaussian 
quadrature, only the order is specified and the points 
and weights are obtained. In 3d, this may be very cum- 
bersome so we start with points symmetrically placed 
and we will show that they give a good approximation. 
Take points on each axis and on the diagonals of the 
cube, (±a, 0, 0), (0, ±a, 0) , (0, 0, ±a) and (±6, ±b, ±b) as 
illustrated in Fig. 13. By symmetry, there are only two 
weights wi and W2 • Thus the integral is approximated by 



1 



dxdydzf{x,y,z) 



{2hf j-hJ-hJ-h 
« wi [f{±a, 0, 0) + /(O, ±a, 0) + /(0, 0, ±a)] 
+ W2/(±6,±6,±6). 



(A2) 



To determine the numbers a, 6, wi and W2 we require 
that every polynomial of order five or less should be inte- 
grated exactly by this scheme. In 3d, this corresponds to 
many different polynomials but we have only four unk- 
nowns and apparently too many equations. This is where 
symmetry comes into play. First, since we integrate over 
a cube from —h to h, all odd polynomials integrate to 
zero. Also, for example, a polynomial of the form x^y^ is 
equivalent to y'^z^, or is equivalent to and and so 
on for every type of polynomials. Thus we only have to 
consider four different polynomials i.e. 1, x^, x'^ and x^y"^. 

Taking f{x,y,z) = 1, the integral gives one and 
thus we obtain the first equation 

6m;i + 8w2 = 1 (A3) 

Taking f{x,y,z) = x^ , the integral gives 

/ dxdydzx^ — ^ and we obtain 



2a^u;i + %h^W2 = — 
o 

Similarly, with f(x,y,z) = 
^2fc)a / dxdydzx'^ = ^ and we find 



(A4) 
we have 



2a^«;i -I- 8b^W2 



Finally, taking f{x,y,z) = x'^y'^, J dxdydzx'^y 



(A5) 

2 _ 



.4 

^ and we obtain 



W2 = 



(A6) 



Solving these equations, we obtain 

40 

Wi = 
W2 = 

a = 
b = 




(A7) 



To make the method adaptive, we take our cube and 
split it in eight. If we take one of these cubes, and put 
the origin in its center we now have the integral over 
a cube from — to ^ centered at Tq. We can thus use 
Eq. (A2) but with /i — > -I and the points {xq ± a, yo, zq), 
{xo, yo±a, Zq), (xo, yo, zo±a) and {xo±b, yo±b, zo±b). We 
can do this for each of the eight cubes obtaining the new 
approximation for the integral I — ^ X)i ^i- "^^^ process 
can be repeated. Each of the eight cubes can be subdivi- 
ded again with integrals from —j to j. When one sub- 
division has converged, this part is stopped. The calcula- 
tion has converged when all subdivisions have converged. 

The convergence criterion requires a detailed discus- 
sion. Assume that we aim at a relative error e. For an 
adaptive method, we need absolute error. Indeed, with a 
simple Id adaptive integration method where one subdi- 
vides the interval in half, if one wishes an absolute error 
5 one usually imposes absolute error | on each of the 
two sub-intervals. To determine the absolute error in our 
case, we first estimate the value of the integral by perfor- 
ming three subdivisions, i.e. using 8^ = 512 cubes. Let 
us call the resulting integral Then we take for the 
absolute error needed for the adaptive integration me- 
thod S = €l^^\ If we imagine launching the integrator 
from scratch, one would ask for ^/8 accuracy in each of 
the 8 sub-cubes when we do a division. This often leads 
to a final answer that is more accurate than desired. For 
heavy numerical calculations it is desirable to optimize 
the choice of the error in each subinterval to minimize 
the computation time while maintaining the final desired 
accuracy. In our case, we claim that it suffices to require 
the absolute error within each subinterval to be ^ ins- 
tead of |, as we might have naively expected. Indeed, if 
we consider each value on the sub-cube as a random va- 
riable li and want an error 6 on the original cube i.e. on 
the sum / = J^i -^i) error for each sub-cube suffices 
is we assume that the errors on the li's are independent 
and uniformly distributed (lUD). Indeed, in that case 
Var(/) = Ei Var(/i) and thus 5i = VSSi,. This hypothe- 
sis of an lUD is of course not rigorous, but we have exten- 
sively tested this choice for the error with many different 
integrands with known integrals. By taking advantage of 
the fact that statistical hypothesis on the errors become 
reasonable since the integral is high dimensional, our ap- 
proach is faster. We also checked that our approach is 
more precise and faster than using three adaptive Id in- 
tegrators. Additional speedup can be obtained by taking 
into account the symmetry of the integrand. 



Annexe B: IPT-D implementation 

In this section we detail how we implemented IPT- 
D. IPT as an approximative solver is fast, but it 
needs to also be implemented in the fastest possible 
way. In IPT, we are solving a system of two nonli- 
near equations with two unknowns : /Uq and /U. The 
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G<,(/-V) = 



j'vbVi + /ill - A(Lij„) 



Repeated until /i"*^' 
= n and D"™=D 



rro = 2Go(r = (r 


.1 »(2-.t) 


(l-f)C/ + /i«-/j 


' no(2 - no) 













E(!W„) = C/- + 



Gf = 



2 1 - BE(2)(!W„) - A(!tJ„) + /.I - S(saJ„) 



updated 



: 2Gf{T = 0-) W'"" = |>^e'-°*E(,u,'„)G/(ia.„ 



calculated numerically and we now explain how to do 
it efBciently. The first necessary step is to calculate 
Go(t). We can show that the function A(ia;„) for a 
one band model with a dispersion relation Sk behaves 
asymptotically like A(iw„)„^oo ^ j^Y^k^l = 
similar to what was previously obtained'^^. For the 3d 
FCC lattice with nearest- neighbor hopping c — 12P. 
We can use this to define G™^ iiuJn) — „-^. , — that 

gives the asymptotic high-frequency behavior of Gq. 
We need it because the Fourier transform necessary to 
get Go(t) must be approximated by a finite sum and 
thus the function must be convergent at least as fast 
as /- ^ yi , while Gn(ia;„) J—. We thus consider the 

function F — Go — Gq"^ instead and add the missing 
terms analytically. 



Figure 14: Flow chart for the impurity solver loop in IPT. 
There is also an outer loop for A, see text. 

equations are n — 2-g ^„ e'""''^G(iaj„) = and D — 

V e'"""^S(ja;„)G(icL;„) = where n and D are fixed 
numbers for a particular set of parameters. The self- 
energy must be calculated using Eqs. (4), (5). We show 
how to do this efficiently. 

We first start with guesses for the hybridization func- 
tion A(ia;n), and the two chemical potentials /^o and ^. 
Then we calculate the impurity model loop as shown 
in Fig. 14. Once the loop has been converged, we use 
the self-energy to calculate the lattice Green's function 
G(iw„) = Y.k »^,.-(e.-'p)-s(»c.„) - With it, the new hybri- 
dization function can be calculated A(ia;„) = — I](ia;„) — 
G~^{iujn) + i^n + A* and finally A, /i and /xq are fed back 
to the impurity model loop. This is repeated until global 
convergence is reached. 

It must be specified here that A is a function of both /i 
and /io but in the present algorithm, once A is fed to the 
impurity loop, it is considered to be independent while 
/i and /io are iterated until we obtain the correct n and 
D. However, once the system is close to convergence, the 
difference between A in the inner loop and the correct 
A(/i,/Lto) becomes really small and once convergence is 
reached, it is indeed the same function. This is also the 
approach that was adopted originally in [12]. This ap- 
proach is much faster than fixing /i and /io, converging 
the entire DMFT calculation, calculating the new n and 
D, updating /x and /io and converging again the DMFT 
calculation until we obtain the correct n and D. We have 
tested and used the two methods and they give the same 
results. The fast method can, for some particular para- 
meter set, become unstable and thus, in these cases, the 
long, more rigorous method, may be used. 

To perform this calculation, we see from Fig. 14 that 
we must calculate Fourier transforms from Matsubara 
frequencies to imaginary time and back. This must be 



N/2-1 

Go{r)=T ^ e-'^(2n+ib7iv^(^^^^ 

-N/2 

-fT^e— "^G™^(zc.„) 

n 

N-l (Bl) 
definition of FFT 

-1-r^e— "-G™^(zc.„). 

n 

In the first sum, imaginary time has been discretized in 
N bins so that j/N — r/fi. There is then a maximum 
and minimum Matsubara frequency. In the second equa- 
lity of Eq. (Bl), FFT stands for Fast Fourier Transform. 
The last term, T^)™ e"*""^Gj,"^(ia;„) can be calculated 
analytically using complex analysis. We find, for t > 0, 

Zi ~ Z2 
Z2 - Zl 

(B2) 

where /(z) is the Fermi function and zj = Mo±-\/Mq +f£^ 
Hence, in terms of FFT's, we obtain 

Go(r) = e— ^"(1/^-1) iFFT(i^(z^„_^./2)) + H>{t). 

(B3) 

We must remember that FFT does not give the value at 
T = /?. We will come back to that point later. We also 
need Go(— t). This can be done using the antiperiodic 
property or by using a procedure similar to Go(t). In 
that case, we obtain (again for t > 0) 

Go(-r) = e— ^(1/^-1) iFFT(i^*(zc.„_^./2)) + H<{t), 

(B4) 
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where 



H<{t) = -^^f(zi)e'^^ + -^^f{z2)e'-^. (B5) 



Zl - Z2 



Z2 - Zl 



With these two Green's function we can calculate the 
second order contribution to the AIM self-energy that ap- 
pears in IPT Eq. (5). In this equation, we need to perform 
an integral. To obtain an asymptotic behavior in Matsu- 
bara frequencies that decays instead of being periodic, we 
cannot do a direct integration^. The trick here is to per- 
form a cubic spline interpolation of 7(t) = Gq(t)Go(— r) 
and then Fourier transform that spline interpolation. To 
obtain the spline, one needs two conditions to solve the 
system of equations. In our case, it suffices to find the 
derivatives at r = 0+ and t = f3~ . We will show later 
that making the Fourier transform of the spline is really 
accurate and introduces a minimum of numerical errors. 

Up to now, we have only considered r > and thus 
what we need are the derivative at 0+ and at . If, for 
the moment, we consider a paramagnetic system we can 
write /(r) = Gq{t)Go{—t). The derivative is thus 

''^^^ -2Go(r)Go(-r)^+G^(r)^^. (B6) 



dr 



dr 



dr 



We calculate the derivative of the Green's functions from 
their definition 

n 
n 

(B7) 

By defining = (-ia;„)F(iu;„) and i^2(«w„) = 

{—icOn)F*{iuJn), we obtain 

= e--^[i/^-iliFFT(fi(ia;„)) 



Zl - Z2 



-f{-zi)e-'^- + 



Z2 - Zl 



-f{-Z2)e- 



^^^izl) ^[1/^-11 iFFT(f2(ia;„)) 

dr p 



+ 



Zl — Z2 



-f{zi)e''^ + 



Z2 — Zl 



-f{Z2V 



(B8) 



With these two equations, we can get the derivatives at 
r = 0+. To obtain the Green's function and its derivative 
at T = I3~ , we use the spectral representation of the 
Green's function to show that Go{f3~) = — 1 — 00(0+), 



dGo(- 



dT 



r=0+ 



and 



T=0+ 



dGo(-T) 
dr 



dr ~ dr 

-/So - 

We now have everything we need to calculate the deri- 
vative of Eq. (B6) for t = 0+ and t = j3~ . Knowing /(t) 

and ^^-^ we can calculate the coefficients of the spline. 



We need the Matsubara-Fourier transform of functions 
represented by that spline on the right-hand side of 



f 



dre" 



7(t 



(B9) 



Let us call S{t) the piecewise cubic spline for /(r). The 
method is presented in details in Appendix E of [32]. 
Integrating by parts, the result for fermionic frequencies 



is 



-Si{Q)-Sn{P) , S'i{Q) + S'^{p) 



+ 



-SnO)-SUfi) 



(BIO) 



(l - e*"™ 
+ 7V-^- -3-^IFFT 



(e^5-r) 



where IFFT is the inverse Fast Fourier transform. This 
result is what we needed since it has the correct high- 
frequency behavior where in principle the first three 
terms are exact while the last one is obtained from a 
numerical inverse Fourier transform. Since the latter is 



the coefficient of 



errors do not adversely affect 



the high-frequency behavior. 



Annexe C: Calculation of iVo(e) and X{e) 

Even though for the DMFT iterations we found 
that the way to obtain accurate results was to use the 
adaptive method described in Appendix A, we show 
here how to calculate No{£), the non- interacting density 
of states and X{s). The latter quantity appears in 
calculations of the conductivity and transport properties 
in general and it is in this context that we used the 
results presented here. Both quantities have the general 
form J2k -F'(k)(5(£ — £/). Since our band structure is not 
simple, we cannot perform the integral analytically. The 
question is thus, how do we treat the delta function in a 
numerical calculation ? 

A simple approach would be to replace the delta 
function by a Lorentzian and perform the integral using 
an adaptive scheme. But, this approximation for the 
delta function gives tails at the edges of the band. A 
Monte Carlo scheme is preferable not only because the 
sharpness of the delta function is maintained without 
tails, but also because one does not need to do a triple 
integral for each value of s : The complete function of e 
is be obtained at once. 

We first choose for how many energy e points we 
want to know the function. This number defines a 
number of bins, shown as dashed lines in Fig. 15. One 
then generates a random point {krc,ky,kz), calculates 
Sk and locates the bin where this number belongs. For 
example, in Fig. 15, the random belongs to the bin n. 
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Figure 15: Illustration of the Monte Carlo integration scheme 
used to obtain No{e) and X{s) 



We then add to this bin ^F(k). This is equivalent to 

approximating the delta function by a rectangle of finite 
width. We continue this process M times and divide, at 
the end, the numbers in the bins by M. The function 
is thus given by the numbers in the bins, each bin 
corresponding to a particular energy e. Accuracy and 
smoothness can be improved by increasing the number 
of random points. 



while for O; = 

t^i^^ ^ (D4) 

+ y (^n - Tn-l) + dn {Tn - . 

with a„, bn, Cn and dn the coefficients of the cubic spline. 

Wo also need to go from Matsubara frequencies to ima- 
ginary time for Gkir) and for F{t) in Eq. (D2). We pro- 
ceed to obtain the analog of Eq. (B3). Ginf has a similar 
structure since the asymptotic behavior of the self-energy 

is Ti{iujn) = Un/2 + ^^—^j— ^- This time the poles are at 

_ (£fc - M + Un/2) ± ^(£fc - + t/n/2)2 + 
2 

(D5) 

where c = f/^n(2 - n)/4. 

Once again, the values at r = /3~ are found using 
the spectral representation of Gk- The expressions are 

GkW-) = -1 - Gk{0+). (D6) 



Annexe D: Optical conductivity 



Gfc(-r) 



= 1 - Gfc(-r) 



r=0+ 



(D7) 



For the optical conductivity, we need the current- 
current correlation function in bosonic Matsubara fre- 
quencies. In DMFT, vertex corrections vanish^, hence we 
have 

Xii(ifi;) = - ^vl^^Ga{k,iuJn)Ga{k,iu>„+ifli), 

(Dl) 

where O; are bosonic Matsubara frequencies while w„ are 
fermionic. To compute this, we again use the convolution 

theorem and FFT, as described for IPT in Appendix B. 
The above equation can be written as 



XiM^i)^- / dTe'''^-y^vlG„{k,T)G. 



.(fc,-r) 



Jo 



dre'^'^Fir). 



(D2) 



The only difference with Eq. (B9) is that here we have 
bosonic frequencies. Using again the cubic spline trick we 
obtain an expression similar to Eq. (BIO). For fli ^ 



+ 



(l - e'"' 



(D3) 



dGkir) 



dr 



dGk{-T) 



{sk-fi + Un/2) ■ 



dGkir) 



dr 



r = + 



(D8) 



dr 



= (£fe - A* + Un/2) - 



dGki-T) 



dr 



=0+ 

(D9) 

We can thus calculate F{t) appearing in Eq. (D2) and 
its first derivative at r = 0+ and t — f3~ to find the cubic 
spline interpolation. 

The sum over wave vectors k requires some comments. 
Like for the DMFT calculation, we use the adaptive 
scheme in Appendix A. However, if we look at F{t) in 
Eq. (D2), in principle for each r we need to perform 
the integral over k independently, but the whole point 
of FFT is to obtain all t points at the same time. Our 
solution is to launch the integrator for all t at the same 
time and keep in memory the estimate of the integral for 
each T while we refine the estimate. Once the integrator 
converges for r = 0, we conclude that this k grid is the 
one for all r and stop the calculation. Since we keep all 
values in memory, we have the function for all r calcu- 
lated for the grid in k space appropriate for r = 0. We 
verified that by converging the calculation for other va- 
lues of T we obtain the same answer. We have also tried 
the other way around, where we interchange the integral 
over k and the Fourier transform in Eq. (D2). In this 
case, we converge the zero frequency and, once again, 
the results are essentially the same. 

Once we have calculated xii(?r2(), we need to ob- 
tain the real frequency representation since the optical 
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conductivity is given by 



x"i (") 



. In the case of CTQMC, 



we use the maxent analytical continuation scheme deve- 



loped in Bergeron et al. [32] for the conductivity. For IPT 
results we use both Pad? and maxent to find xii(w). 
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